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Large eddy simulation of a forced round 
turbulent buoyant plume in neutral surroundings 

By A. J. Basu AND N. N. Mansour 
1. Motivation and objectives 

Buoyant flows play an important role in various technological and environmental 
issues. For example, dispersal of pollutants, smoke, or volcano exhaust, in the 
atmosphere, vertical motion of air, formation of clouds and other weather systems, 
and flows in cooling towers and fires are all determined primarily by buoyancy 
effects. The buoyancy force in such flows can originate from either a heat source or 
due to different densities between a fluid and its surroundings. Whatever the cause, 
the flow can be understood by studying the effects of the tight coupling between 
the thermal and the velocity fields since density differences can be characterized as 
temperature differences. 

Heat addition to turbulent flows always leads to very complicated behavior, 
greatly affecting their structure and entrainment characteristics. Heat, can be added 
in many different ways into a flow; for example: at the source, or off-source. These 
two methods of heat addition can have very different effects on the development 
of the flow. Basu & Narasiinha (1999) studied the effects of off-source volumetric 
heating (similar to that due to latent heat release in a. cloud, for example) using 
Direct Numerical Simulation (DNS) of a. circular jet-like flow and found that the 
large-scale structures break down and entrainment is inhibited. Here we look at a 
round turbulent plume with both momentum and buoyancy added at the source, 
using the tools of Large Eddy Simulation (LES). 

Similarity solutions of turbulent circular plumes have been proposed in the past 
by Zeldovich (1937) and Batchelor (1954) based on the Boussinesq hypothesis, which 
assumes that the density changes in the flow are small compared to the ambient 
density, an assumption that is valid in regions away from the plume source. In a 
neutral environment (no change in ambient density with height), the mean vertical 
velocity and buoyancy acceleration are given respectively by: 
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where W is the mean vertical velocity along the axis of the plume, F 0 is the rate 
of addition of buoyancy, c is the height (from the buoyancy source), g is acceler- 
ation due to gravity, Ap is the density difference at a point with respect to the 
ambient density p a , and A\v,B\\' are the parameters which quantify the Gaussian 
fit to the mean velocity profile, while Ar.Bj ■ arc' the corresponding ones for the 
density (or, equivalently, temperature) profile; // represents the similarity variable 
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i] = r/z, where r is the radial distance at any z. For a flow involving sources of 
both momentum and buoyancy such as the one studied here, it is known that the 
buoyancy effects overwhelm the momentum added at the source, typically after an 
axial distance / L ; \ / > 5. where 

Lm = M^/Fy* ( 3 ) 

is the so-called Morton lengthscale. The source momentum flux M 0 and buoyancy 
flux F 0 are defined as: 

fRo 

M 0 = 2 tt / W 2 rdi\ ( 4 ) 
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F 0 = 2n / Wg rdi\ (5) 
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where R 0 is the source radius. 

The forced buoyant plume has been the object of many experimental investi- 
gations in the past (see List 1982, Rodi 198C for reviews). Most notable among 
the reported experiments are the ones lw Rouse et al (1952), Nakagome & Hi- 
rata (1977). Papanicolaou List (1988), and Shabbir & George (1994) in that 
they report widely different results. The primary disagreements among the vari- 
ous experimental results are regarding the centerline values of mean velocity and 
buoyancy profiles as well as the plume spreading rate. For example, the spread 
of reported experimental values (as compiled in Shabbir & George) are given here 
inside brackets for the different parameters for Gaussian fit of profiles described 

above: A w = (3.4 - 4.7), B w = (55 - 96), Ar = (9-1 - 14.28), B r - (48 - 71). 
Clearly the scatter is quite large. 

The discrepancies between different experimental data are likely to be due to var- 
ious factors, some of which are detailed below: (i) boundary effects of the solid wall 
lateral boundaries and presence of reverse or co-flow can influence the entrainment, 
(ii) measurements may not have been carried out in fully developed turbulent region 
in some cases, (iii) hot-wire anemometer measurements are known to be insensitive 
to direction, and therefore the measurements made outside the half-width of the 
plume may not be reliable in such cases. Since F 0 is used to scale for self-similarity, 
a lot crucially depends on the accurate determination of F a , the measurement of 
which, unfortunately, can be influenced by the experimental errors. As will be seen 
later, the above mentioned problems with the available experimental data somewhat 
constrain the validation of numerical predictions made in this study. 

2. Accomplishments 

2 A Governing equations 

In the present study, we aim to compute the evolution of a circular plume in neu- 
tral surroundings with momentum and buoyancy added only at the source. Since 
the density differences away from the source of a plume are known to be small 
compared to ambient density, the Boussinesq approximation is assumed to be valid. 
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whereby the effects of density variations are modeled by a source term in the mo- 
mentum equations. The resulting equations in the non-dimensional form are given 
by: 

V*u = 0, (6) 

| + V.(uu) = -V, + ^u + |lef, (7) 

f + (u v.e^v-’e, (8) 

where u is the velocity vector, p is pressure, 0 — (T — T a )/(T Q — T a ) is the tem- 
perature difference (where T, T a ,T 0 are the absolute temperatures, respectively at 
a given point, the ambient and the plume source 1 ). The non-dimensional governing 
parameters are: 

Re = U 0 D 0 /u = Reynolds number, 

Pr = v/k — Prandtl number, 

Gr = (\(j (T 0 — T a )D'tlv 2 — Grashoff number, 
where U 0 and D 0 are respectively the velocity at and diameter of the source, v 
is the kinematic viscosity of the fluid, k the thermal diffusivity, o the coefficient 
of thermal expansion, and g the acceleration due to gravity. The last term in the 
momentum equation (7) above is the buoyancy source term and acts only along the 
vertical z direction. 

2.2 Numerical technique 

Following the standard approach for LES, the above equations are passed through 
a spatial filter, and the resulting ('fleets of the subgrid scales are modeled. In 
the present study, we have used the well-known dynamic model for this purpose 
(Germano et ai 1990, Lilly 1992). 

The spherical polar coordinate system ( r, 0, 0, along the radial vector, lateral, and 
azimuthal directions respectively) is used here since, for the present flow with its 
conical mean growth, a spherical system allows for a well-balanced resolution of the 
flow field with a reasonable number of grid points. The present computer code is a 
modified version of that developed by Boersma et al. (1998) for cold jets. The basic 
numerical formulation being similar, the reader is referred to the above-mentioned 
paper for various implementation details. We will only describe some fundamental 
features of the numerical scheme and also the extensions made in the present study. 

Second order accurate finite volume technique is used to approximate the equa- 
tions of motion in space, along with a second order explicit Adams-Bashforth scheme 
for time-integration. Pressure correction technique is employed after each Adams- 
Bashforth step in order to ensure that the velocity field remains divergence-free at 
all times. The resultant Poisson equation is solved using fast. Fourier transform 
along the periodic (0) direction and cyclic reduction along the other two directions. 
The convective term in the temperature equation is treated using a TVD scheme 
(see, e.g., Vreugdenhill &: Korcn 1993) in order to keep 0 between the bounds 0 
and 1. 
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2.3 Boundary condition* 

The flow considered, as mentioned earlier, is a forced circular buoyant plume 
issuing from an orifice in a wall at the bottom end of the computational domain; 
we assume a topliat velocity profile. The boundary conditions at the inflow are 
therefore: 

u r — W 0 cos6 in the orifice and u r — 0 elsewhere, 
uq — — W 0 sind in the orifice and a# = 0 elsewhere, 
e (j> — 0? 

0 = 1 in the orifice and 0 = 0 elsewhere, 

where W 0 = 1 is the non-dimensional axial velocity in the orifice (diameter D 0 = 1). 
Note that u r , uq, here represent velocity components along the r, 8, (f> directions 
respectively. 

At the lateral boundary, the so-called traction-free boundary conditions are used 
(Gresho 1991 ): 

(?ij • 7i j = 0 , ( 9 ) 

where a, 7 is the stress tensor, and n, the unit normal to the boundarv. The ad- 
i j j 

vantage of this boundary condition is that entrainment is allowed (Boersma et al. 
1998). For temperature, the normal derivative is set to zero at the lateral boundary: 


oe 

On 


= 0. 


( 10 ) 


At the outflow boundary, we use the so-called advective boundary condition, as 
in Boersma et al. . Thus, for any velocity component u, we have 


Dt A dr' 


(ID 


where the advective velocity W,\ is a function of 6 and is obtained at each step 
by averaging the streamwise velocity along <p near the outflow domain and setting 
negative values of Wa to zero. As mentioned in Boersma et al. , this outflow 
boundary condition is not entirely satisfactory for the present flow. Therefore, 
we also add a so-called “buffer zone 11 near the outflow region where terms similar 
to — A(u — u t. ar gct) are added to the momentum equation in order to damp the 
turbulence in this region; here A is a space- varying positive parameter (A = 1 at the 
outflow boundary and decays to 0 outside the buffer zone), and u target a desired 
velocity field. Typically, an analytically or experimentally derived mean velocity 
field near the outflow domain is used as u target- 

2.4 RrsulU 

The parameters defining the flow in the present computation are chosen so as 
to closely correspond to the experimental situation for the round turbulent buoy- 
ant plume reported in Shabbir Sz George (1994). Therefore, we use the following 
non-dimensional numbers that completely determine the flow: Re = 3500, Pr = 
0.7, Gr = 8.575 x 10 b . White noise with a peak amplitude of 0.02 has been added at 
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all times at the source in order to mimic the disturbances in the exit profile reported 
in the experiment. 

The computations are carried out in a domain that extends to 50 diameters down- 
stream of the source. The domain encompasses a conical volume of lateral angle 
7 t/ 15, with a virtual origin that is 15 diameters upstream of the orifice. This volume 
is discretized using a grid of size (N r = 128. Ng — 40, A r ^ = 32). where N r , Ng, 
are the number of finite- volume cells along the (r, #,<£>) directions respectively. The 
grid spacing along r increases linearly from 0.1122 near the source to 0.069 near the 
outflow boundary. Grid spacing is constant along 8 and <6. 

The computations are carried on for 70,000 time-steps in the present study. The 
time taken on a 8-processor SGI Origin 2000 is of the order of 1 second per step. 

In the results presented below, the data from the last 10 diameters are not shown 
since the effects of outflow boundary conditions and buffer zone are likely to be 
significant there. 

2.4 -1 Temperature distribution 

The distribution of non-dirnensional temperature difference 0 in a buoyant plume 
is an important parameter since it represents the driving force behind the* flow, espe- 
cially after a few diameters downstream of the orifice where' the effect of initial mo- 
mentum becomes insignificant in comparison with the buoyancy force. I 11 Fig. 1(a), 
we present a. contour plot of the typical instantaneous distribution of 0 over a ver- 
tical section passing through the axis of the plume. The contour levels go from 0.1 
to 1 in steps of 0.1. Very clearly, high values of 0 are limited to the region close 
to the orifice. Temperature' differences higher than 0.2 (the non-dimensionalized 
temperature difference at the source being 0 = 1) exist only within the first 10 
diameters or so. After about 20 diameters from the source, peak 0 anywhere in the 
plume has fallen to 0.1 or below, indicating that the Boussinesq approximation is 
valid over the self-similar region of the flow, the region of our primary interest. 

The lower levels of 0 (those between 0.01 and 0.1) at the same time are high- 
lighted in Fig. 1(b) in order to complete' this instantaneous picture of temperature 
distribution. The important point to note is that there does not appear to be any 
“holes 1 " in the temperature distribution, unlike' those reported for passive scalar 
(Papantoniou & List 1989). Visual analysis of data indicated similar qualitative* 
nature of distribution at other times. 

2.4-2 Streamwise variation of momentum and buoyancy fluxes 

The similarity solution for plumes (e.g., Batchelor 1954) suggests that, for an ideal 
gas plume in a neutral environment, the buoyancy flux is conserved. The momentum 
flux, in comparison, grows with streamwise distance due to entrainment. Figures 
2(a), for momentum flux, and 2(b), for buoyancy flux, show that this predicted 
behavior is observed in the presently computed results. Note that the quantity F 
in the figure* represents 


F = 2n 


Ho 


Jo 


(\g(W AT+ < irt > )rdr. 
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FIGURE 1. Distribution of computed temperature difference 0 in the buoyant 
plume at a non-dimensional time of 150: (a.) contour levels from 0.1 to 1 in steps of 
0.1, (b) contour levels from 0.01 to 0.1 in steps of 0.01. The computational domain 
boundaries along the lateral direction are also shown. 

where < irt > represents the turbulent contribution to the total buoyancy flux, te 
and / being the fluctuations in the streamwise velocity and temperature respectively. 

The buoyancy flux shows a fall of about 15% in its value (with respect to that 
at the source) over the first 5 diameters as we go downstream from the source but 
remains quite constant afterwards. This indicates that the present computation 
quite faithfully represents the integral behavior of a buoyant plume in the self- 
similar region downstream, where the Boussinesq approximation is valid. The order 
of change in F over c predicted here is similar to that reported in the various 
experiments. 

2-4-3 Mean flow 

The results for streamwise variation of centerline mean velocity W r and tempera- 
ture difference A T r are shown using the usual non-dimensional form (see, e.g., Shab- 
bir George 1994) in Fig. 3 using a log-log plot. As can be seen very clearly, the 
temperature distribution takes longer to achieve self-similarity when compared to 
the streamwise' velocity; thus W (: achieves self-similarity around z/L\j — 6, whereas 



LES of a buoyant plume 


7 




FIGURE 2. Streamwise variation of radially and azimuthally integrated (a) mo- 
mentum and (b) buoyancy fluxes. 

AT r does so by zjL\\ = 15. The scatter in the various reported experimental data 
is quite large, and the present results fall towards the high end of this scatter. The 
computed streamwise growth rates, however, are remarkably consistent with that 
predicted by the theory and experiments. 

Next we look at the mean radial profiles of < W > and < AT > at some 
streamwise stations in order to establish the self-similarity of the present computed 
results. Figure 4 shows the computed results along with Gaussian fits of available 
experimental data as compiled by Shabbir k George. The computed velocity profile 
closely matches that given by Rouse et al. (1952). is similar in width compared to 
the ones in Papanicolaou k List (19SS). but is different from those reported by 
Shabbir k George (1994) and Xakagome k Hirata (1977) both in terms of the 
width of the profile and the maximum value at. the centerline. 

The computed temperature profiles differ from all the experimental ones in that, 
the present ones are narrower. The centerline values, however, are close to that 
reported by Papanicolaou k List. 

2-4-4 Turbulence properties 

In Fig. 5. we present the computed results for the second order moments of 
fluctuating quantities using solid curves, along with the experimental data from 
Shabbir k George (1994) using symbols. The radial profiles of turbulence stresses 
at different streamwise stations converge quite well in the time over which they 
are averaged and show reasonable collapse when plotted using similarity variables. 
The present results for < w 2 > and < uw > are comparable to those obtained by 
Shabbir k George except for the width of the profiles. For the fluctuating stresses 
involving temperature, however, it is the magnitude which shows large differences, 
consistent with similar relative differences seen for radial mean profiles in Fig. 4 
above. However, a few words of caution are in order: such comparisons between 
filtered stresses (as in the present case) with experimental results may not be very 
meaningful unless it can be shown that the subgrid stresses are small in comparison. 
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FIGURE 3. Normalized mean (a) buoyancy and (b) vertical velocity along the 
plume centerline plotted against the normalized distance z/L\j from the source. 
The bold lines represent the present computation. Fits to various experimental 
results are shown using lines with symbols (□ : Shabbir & George 1994, A : Papan- 
icolaou L List 1988, v : Rouse at al 1952, o : Nakagome & Hirata 1977). 




FIGURE 4. Radial profiles of (a) mean streamwise velocity and (b) mean temper- 
ature difference plotted in similarity form. The computed results at zjD — 15 to 
40 at interval of 5 are shown using bold lines. The' line's with symbols represent, the 
experimental data as described in Fig. 3. 

3. Conclusions and future plans 

In the present preliminary study, we have' shown that reasonable predictions of 
the evolution of a turbulent round buoyant plume can be obtained by means of LES 
with the dynamic subgrid model. The large scatter between available experimental 
data makes definitive comparison difficult. However, the present results fall within 
the scatter region of the experimental data. 

Judging from the scatter in various experimental results, the boundary condi- 
tions are likely to play a significant role in this flow. The present computations 
also bring up various issues regarding the imposed conditions along the lateral and 
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FIGURE 5. Second order moments plotted in terms of similarity variables, (a) 
the vertical velocity fluctuations, (b) the turbulent shear stress, (c) temperature 
fluctuations, (cl) the vertical turbulent heat flux. The lint's represent the computed 
results at zjD - 20, 25, 30, 35 and 40: the experimental data from Shabbir & George 
(1994) are shown superposed using □ symbols. 


outflow boundaries. The traction-free condition along the lateral boundary can in- 
fluence entrainment, and so can the advective condition and buffer zone along the 
outflow boundary. What should be the nature of the physically correct boundary 
conditions along such boundaries still remains an open question (see Gresho 1991 
in this regard). We plan to probe further into this aspect, and study the effect of 
different boundary conditions. The effect of grid resolution and position of lateral 
boundaries will also be investigated. 

Long term plans include LES of plumes in stratified environments and inclusion 
of cloud models in such flows. 
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